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ABSTRACT 



Following Papadakis (2005) 's numerical exploration of the Chermnykh's 
problem, we here study a Chermnykh-like problem motivated by the 
astrophysical applications. We find that both the equilibrium points and 
solution curves become quite different from the ones of the classical planar 
restricted three-body problem. In addition to the usual Lagrangian points, 
there are new equilibrium points in our system. We also calculate the Lyapunov 
Exponents for some example orbits. We conclude that it seems there are more 
chaotic orbits for the system when there is a belt to interact with. 
Subject headings: planetary systems - stellar dynamics 
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1. Introduction 

Chermnykh's problem concerns the motion of a test particle in the orbital plane of a 
dumb-bell. This dumb-bell rotates with constant angular velocity n around its mass center. 
The equations of motion of this test particle are (Chermnykh 1987): 

dx 

— — u 
dt 

dy 

— — v 
dt 

du _ ^ dU* 

dt dx 

dv n dU* 

— = —2nu 



dt dy 



where the potential U* is 



^ = -V^ 2 + ^----' ( 2 ) 
2 n r 2 



T\ = \J(x + /i 2 ) 2 + y 2 , r 2 = \J (x — Hi) 2 + y 2 and /ii + ii 2 — 1- Usually, the mass parameter 
/i = /i 2 and the parameter n G [0, oo). 

This interesting problem was first analysed in Chermnykh (1987) on the Lyapunov 
stability of the triangular solutions. Gozdziewski and Maciejewski (1998) studied the 
Lyapunov stability of the Lagrangian points in the entire range of angular velocity n and 
mass parameter \i. Recently, Papadakis (2005) investigated the equilibrium points and 
the zero-velocity curves of this problem. In particular, some families of periodic orbits are 
studied in great detail. The stability zones are also discussed in that paper. 

In general, this important problem has many applications in celestial mechanics (see 
Gozdziewski and Maciejewski 1999) and chemistry (see Strand and Reinhardt 1979). . 
When the angular velocity n — 1, it is the classical planar restricted three-body problem. 
When n = 0, we get the Euler's problem of two fixed gravitational centers. 
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On the other hand, because (i) there are many discovered extra-solar planetary systems 
(Laughlin & Adams 1999, Rivera & Lissauer 2000, Jiang & Ip 2001, Ji et al. 2002), in 
which about five of them are in binary stars and, (ii) the effect of belts might be important 
according to some of our previous studies (Yeh & Jiang 2001, Jiang & Yeh 2003, 2004), we 
here construct a Chermnykh-like problem motivated by these astrophysical applications. In 
principle, we consider the angular velocity n > 1 and its value is determined by the belt's 
gravitational potential, which is considered as part of the total potential to govern the test 
particle's motion. 

Therefore, we construct our basic model in Section 2. In Section 3, we study the 
existence of new equilibrium points. We discuss the behavior of the orbits with initial 
conditions near the new equilibrium points in Section 4 and calculate their values of 
Lyapunov Exponent in Section 5. Section 6 concludes the paper. 



We consider the motion of a test particle influenced by the gravitational force from the 
central binary and the circumbinary disc. 

For convenience, we set the gravitational constant G — 1. We assume that two masses 
of the central binary are /ii and /i 2 and choose the unit of mass to make /ii + /i 2 = 1. The 
separation between two stars in the central binary is set to be unity. The time unit is 
therefore determined from the above choice. 

Under the influence from the belt, when ^ — /i 2 — 0.5, both components of the 
central binary move on circular orbits at r = 0.5 with mean motion, i.e. angular velocity 



2. 



The Model 



n = 




1 — 2/ 6 (0.5). The gravitational force fb from the belt is 
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where 



= 

_l_ ty> 



F{£) and are elliptic integral of the first kind and the second kind and 

rn/2 I 

F(0 = / 

Jo ^/l - £ 2 sin 2 0' 

r7 r/2 



E(0 = f ^l-^sinV^'- (5) 

(Please see Appendix A for the derivation of Eq.(3).) Hence, the x and y components of 
gravitational force are 

_dV_ _ x 

o Jb j 

ox r 

Jb j 

dy r 

(6) 

where /& is in Eq.(3) and F is the potential from the belt. 

Because both components of central binary are doing circular motions, the equation 
of motion can be written on the frame rotating with the central binary. The equations of 
motion of this problem are 

dx 

— — u 
dt 

dy 

— — v 
dt 

du_ 2nv _ dlP _ dV 
dt dx dx 

dv _ dU* _ dV 

dt ~ dy dy ' 

(7) 

where the potential U* is given by Eq.(2). r\ = \J(x + /i 2 ) 2 + y 2 and r 2 = <J(x — /ii) 2 + y 2 . 
The Jacobi integral for this system is therefore 

Cj = -u 2 - v 2 - 2U* - 2V. (8) 
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The density profile of the belt is 





p(r) 



when r < r\ or r > r Q , 



c 



c 



{ C0S [ffc^)] } WlleI1 r i< r < r a, 

when r a <r < r b , 
k ^Hlffeft]}' whenr 6 <r<r , 
where r = y/x' 2 + y' 2 , c is a constant completely determined by the total mass of the belt 
and p is a natural number. In this paper, we set p = 2, r a = r» + 0.1, rj, = + 0.9, and 
r o = r i + 1-0 for all numerical results (Lizano & Shu 1989). It is a power-law profile with 
smooth edges. Hence, the total mass of the belt is 

r2ir rr, 



(9) 



M b = ^ J p(r')r'dr'd<p 



= 2txc 

To 1 



+ 



rr„ I 



r a 1 
COS 



COS 



7r (r 



2 (r, - r a ) 



dr' + ln(r 6 /r ) 



7r (r - r b ) 



dr' 



(10) 



.2 (r -r 6 )_ 

We use r\ and M b as the controlling parameters of the belt profile. Using Eq.(8), we 
determine the zero-velocity curves of our system when Ti = 0.7 and M b = 0.3, as shown 
in Figure 1. Because there are elliptic integrals in the potential V, this has to be done 
numerically and the points are not uniformly distributed on the curves. For the results in 
Section 4 and Section 5, we always set = 0.7 but M b = 0.3 or 0. Following Tancredi et al. 
(2001), we will integrate the orbits up to 1000 binary period. 



3. The Equilibrium Points 

We know that there are five equilibrium points (Lagrange points) for the classical 
restricted three-body problem. Thus, our system shall have five equilibrium points when 
M b = 0. It would be interesting to investigate the number of equilibrium points when 



-7- 



Mb > 0. When there are more than five equilibrium points, we claim that the new 
equilibrium points exist. 



The equilibrium points (x e ,y e ) of System (7) satisfy the following equations: 
f(x,y) = n 2 x- 

g(x,y) 



n 



2 y\i\ y^2 . y , , , n 
n v- — - 73- + -A ( r ) = o, 
"2 ' 



Because we consider ^ — /i 2 — 0.5 here, we have: 

(a; + 0.5) (a; -0.5) x 



n 2 x — 



f(x,y) 
g(x,y) = n 2 y 



2rf 



Irk r 



y y . y , , , n 



+ -/ 6 (r) = 0, 



where T\ — \J (x + 0.5) 2 + y 2 and r 2 = \J (x — 0.5) 2 + y 2 . 
For convenience, we define 



Mi/) = ^ - 



- 3/2 + AW 



(0,2/) 



and 



(x,0) 



k(x) = n z x ^ — ^ - + -f b (r) 

v ; 2|a; + 0.5| 3 2\x - 0.5| 3 r v J 

We have the following properties for the case of equal-mass binary: 
Property 3.1 

If(x e ,y e ) is an equilibrium point of System (7), then we have: 
either (1) x e = and y e satisfies h(y) = 
or (2) x e satisfies k(x) =0 and y e = 0. 

Proof: Suppose that (x e ,y e ) is an equilibrium point, thus it satisfies f(x,y) 
g(x,y) = 0. From Eq.(14), we have 

1 



(11) 
(12) 



(13) 
(14) 



(15) 



(16) 



and 



JJe 



n 



2r? 2r 3 



1 , Mr) 

2 r 



= 0. 



(Ze,3/e) 
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Hence, y e = or [n 2 - ± - A + 



fb(r)] 



(Xe,Ve) 

0: 



0. We now discuss these two cases separately. 



Since (x e ,y e ) is an equilibrium point, that is f(x e ,y e ) = 0, we have: 



= f(x e ,y e )=x e 



n 



2rf 2r? 



1 + Mr) 
■i r 



1 



1 



4rf 4r| 



1 1 

^3 ' 



(17) 



Thus, ri = r 2 , i.e. (x e + 0.5) 2 + y\ = (x e — 0.5) 2 + y\. Hence x e = 0. We have r\ = r 2 



Vl/4 + l/e 2 and thus, h(y e ) = n 2 - [\ + y 2 ]" 3 / 2 + 



MO 



n 



(0,2/e) 



i I Mil l 

2rf 2rf ~ r J 



= 0. 



(0,y e ) 



Therefore, x e = and y e satisfies = for the equilibrium point (x e: y e ). 



(II) Ve = 0: 

f(x e ,y e ) = f(x e ,0) = k(x e ) = for the equilibrium point (x e ,y e ). Thus, x e satisfies 
k(x) = and y e = 0. <0> 

Property 3.2 

(^4^ If y e satisfies h(y) = 0, then (0,y e ) is an equilibrium point of System (7). 
(B) If x e satisfies k(x) = 0, then (x e , 0) is an equilibrium point of System (7). 

Proof of (A): Suppose that y e is one of the roots of h(y) , i.e. h(y e ) = and we set x e = 0. 
Because g(0,y e ) = y e h(Ve) = 0. By Eq.(13), /(0,y e ) = -\[\ + y 2 ]" 3 / 2 + \[\ + y 2 ]" 3 / 2 = 0. 
Thus, (0,y e ) is the equilibrium point of System (7). 

Proof of (B): Suppose that x e is one of the roots of k(x), i.e. k(x e ) = and we set 
y e = 0. Since y e = 0, it is trivial that g(x e , 0) = 0. Because k(x e ) = 0, f(x e , 0) = k(x e ) = 0. 
Thus, (x e ,0) is the equilibrium point of System (7). <0 

Because of the above two properties, in stead of searching the roots for two-variable 
functions f(x,y) and g(x,y) to determine the equilibrium points on the x — y plane, we 
only need to find the roots of one variable functions h(y) and k(x) to get all the equilibrium 
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points of System (7) when the mass parameter /j, = 0.5. Therefore, we numerically solve 
both h(y) = and k(x) = and find out the number of equilibrium points for different given 
parameters. The numerical scheme of root finding is the Van Wijngaarden-Dekker-Brent 
Method (Brent 1973). This is an excellent algorithm recommended by Press et al. (1992). 
We set a high level of accuracy that the maximum error is 1CT 8 for the locations of 
equilibrium points on both x-axis and y-axis. 

Figure 2 shows the results on the r\ — M b plane. That is, we numerically search 
the solution of k(x) = and h(y) = for different (r,,M 6 ), where 0.5 < r\ < 1 and 
0.001 < M b < 0.35. The grid size is 0.01 for and is 0.02 for M b . Those (rj, M b ) with new 
equilibrium points are marked by full triangle points. 

We find that the region with new equilibrium points is itself a triangle. To have new 
equilibrium points, the mass of the belt has to be larger than 0.15 and is somewhere 
between 0.7 and 0.8. 

It is not surprising that there are more equilibrium points in our system because the 
potential field of whole system is more complicated than the classical restricted three-body 
problems and thus there are more possibilities that gravitational forces from different 
components can balance each other. 

4. The Orbits 

Since we have discovered new equilibrium points for particular regions of r\ — M b plane 
for our system, it would be interesting to investigate the orbital behavior around these new 
equilibrium points. 

From Figure 2, we notice that when M b = 0.3, r\ = 0.7, we could have new equilibrium 
points. We thus use this as a standard case. At the y-axis, we explicitly determine the 
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locations of all equilibrium points, which are at (x e ,y e ) = (0, ±0.73), (0, ±0.77), (0, ±1.05), 
in addition to the usual LI point at (0,0). (We find that all new equilibrium points are at 
y-axis because of the choice /ii — /i 2 — 0.5.) 

To understand the properties of these equilibrium points at the y-axis, we numerically 
calculate the orbits with initial conditions very close to these points. There are four initial 
conditions in this paper, i.e. 

(a) Initial Condition LI: (x,y,u,v) = (e, 0,0,0), 

(b) Initial Condition El: (x,y,u,v) = (e, 0.73, 0, 0), 

(c) Initial Condition E2: (x,y,u,v) = (e, 0.77, 0, 0), 

(d) Initial Condition E3: (x,y,u,v) = (e, 1.05, 0, 0), 
where e = 0.001. 

Figure 3 is the orbits on the x — y plane for these four initial conditions with M& = 0.3. 
There are four panels, i.e. panel (a) is for Initial Condition LI, panel (b) is for Initial 
Condition El, panel (c) is for Initial Condition E2 and panel (d) is for Initial Condition 
E3. To have the idea about the disc's influence, we also integrate the orbits with these four 
initial conditions while the disc mass M b = as shown in Figure 4. 

For Initial Condition LI, Figure 3(a) shows that the orbit fills the Roche lobe, i.e. the 
zero-velocity curve passing the LI point. The test particle is gravitational bound to one of 
the star because it probably has less initial total energy than those ones in Figure 3(b)-(d). 
Figure 4(a) shows that the result is similar when M b = 0. 

For Initial Condition El, Figure 3(b) shows that the test particle moves randomly at 
the beginning and then steadily move out of the region of central binary. Figure 4(b) shows 
that the test particle drifts away from the initial location regularly. 

For Initial Condition E2, Figure 3(c) is similar to Figure 3(b) but the test particles 
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move more chaotically in the central region and then move out slowly at the first stage and 
more quickly finally. Figure 4(c) is very similar to Figure 4(b). 

For Initial Condition E3, Figure 3(d) shows that the test particle is always close to the 
equilibrium point and moving on a regular periodical orbit. The equilibrium point could be 
a neutral stable point. Figure 4(d) shows that the test particle is on a complicated orbit 
and drifts away from the central region slowly. 

We notice that the orbits in Figure 3(a) and Figure 4(a) are surrounding one star only. 
Most other orbits show that the test particles move spirally outward from the central region 
except the one in Figure 3(d). 

5. Lyapunov Exponent 

In addition to analyze the orbits themselves, we calculate these orbits' Lyapunov 
Exponents to understand how sensitively dependent on the initial conditions. In general, 
the larger value of Lyapunov Exponent means more sensitively dependent on the initial 
conditions. 

Our procedure to calculate the Lyapunov Exponent was from Wolf et al. (1985). The 
maximum Lyapunov Exponent is defined by 7 = lim^ 00 x(t), where x(t) is the Lyapunov 
Exponent Indicator. The system is chaotic if 7 > or otherwise regular. However, it is 
not possible to take t — > 00 in practice. We numerically determine the Lyapunov Exponent 
Indicator x(t) and usually plot ln%(t) as function of lnt after the evolution is followed up 
to some time. If the curve (or the envelope of the curve) on the lnt — lnx(t) plane shows a 
negative constant slope, the system is regular, otherwise it is chaotic. 

Figure 5(a)-(d) are the results of Lyapunov Exponent for Initial Condition LI, El, E2 
and E3 respectively. There are two curves in each panel, where solid curve is the result of 
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M b = 0.3 and dotted curve is the result of M b = 0. 

From these results, we notice that the dotted curves in Figure 5(b)-(c) and the solid 
curve in Figure 5(d) show negative constant slopes. Thus, when M b = 0, the orbits of Initial 
Condition El and E2 are likely to be regular. When M b = 0.3, only the orbit of Initial 
Condition E3 is likely to be regular. This is probably due to that the equilibrium point 
close to Initial Condition E3 is neutral stable. Note that the solid curves in Figure 5(b)-(c) 
are shorter because the test particles were ejected and the calculations were stopped when 
the distances were larger than 10 4 from the central region. 

We can also notice that, in Figure 5(a), the solid curve is always higher than the 
dotted curve. That is, the values of Lyapunov Exponent Indicator of the orbit with Initial 
Condition LI when M b = 0.3 are always larger than the one of the orbit with Initial 
Condition LI when M b = 0. To summarize, the existence of belt seems to make the system 
to have more chaotic orbits. 

6. Concluding Remarks 

We have provided the equations for a model which includes the gravitational influence 
from a belt around the central binary. We find that, in addition to the usual Lagrange 
points, there are new equilibrium points on the y-axis when each mass of the components 
of the binary is equal to 0.5. To study the orbits around these new equilibrium points, we 
calculate the orbits and their Lyapunov Exponents with four different initial conditions. It 
seems that the system could have more chaotic orbits when there is a belt. 
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A. Gravitational Force from the Belt fb 



Let F(£) be the elliptic integral of the first kind and assume 



( = (AD 
r + r 



By sin 2 A = |(1 — cos 2A), we have 



F(C) = / 

^1-esinV 



o 



(r + r') / . 

JO ^ r 2 + r /2 + 2rr ' C OS(20') 



Jo 



2 7o ^ r 2 + r i2 + 2rr ' COS (e) 

(r + r') f n 1 



2 Jo yV 2 + r' 2 — 2rr'cos(0) 
Therefore, from Equation (A2), we have 



/ , # 
Jo J r 2 _|_ r /2 _ 2rr'cos(0) 

2 / # 

Jo Jr 2 + r' 2 — 2rr J cos(0) 



4F(0 
r + r' 

The gravitational potential is 



V{r) = - ( r pl{r ^ r ' dr'd<j>' 

Jr'JO yV 2 + r /2 _ 2rr / C OS(0') 

= _ 4 , ^(OpV)^ 



(A2) 



(A3) 



Jr' r + r' v ' 



r + r 

where Equation (A3) has been used. Now we differentiate this potential with respect to r 
to get the gravitational force, so 
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= 4/ 

Jr> r + r' 
J r > or \r + r 



dr'. 



(A5) 



Let E{£) be the elliptic integral of the second kind and £' = VI — C 2 - Since is the 
elliptic integral of the first kind, from Byrd & Friedman (1971), we have 



E 



/2 



By Equation (4) and (A6), we calculate 



d£ dr 



E 



r (r + r'y 



£(r + r') F(r - r') 



2 r ( r _ r ') 2r(r + r') 
We substitute Equation (A7) into Equation (A5), so we have 

h = 



dV_ 

dr 



I r ' 



p{r')r' 



rp rp< ty _l_ 



dr'. 



(A6) 



(A7) 



(A8) 
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Fig. 2. — The results of the numerical survey on equilibrium points. Those (r*j, M b ) with 
new equilibrium points are marked by full triangle points. 
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Fig. 3. — The orbits when the disc mass M& = 0.3. Each panel is for different initial 
conditions, i.e. (a) Initial Condition LI; (b) Initial Condition El; (c) Initial Condition E2; 
(d) Initial Condition E3. 
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Fig. 4. — The orbits when the disc mass M b = 0. Each panel is for different initial conditions, 
i.e. (a) Initial Condition LI; (b) Initial Condition El; (c) Initial Condition E2; (d) Initial 
Condition E3. 
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Fig. 5. — Lyapunov Exponent for different disc mass and initial condition, i.e. Panel (a) is for 
Initial Condition LI; Panel (b) is for Initial Condition El; Panel (c) is for Initial Condition 
E2; Panel (d) is for Initial Condition E3. In each panel, solid curve is for M b = 0.3, dotted 
curve is for M& = 0. 



